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Abstract 

We investigate the dynamical evolution of the Orion Nebula Cluster (ONC) by means of direct 
iV-body integrations. A large fraction of residual gas was probably expelled when the ONC formed, 
so we assume that the ONC was much more compact when it for med compared to i ts cur rent size, 
in agreement with the embedded cluster radius-mass relation from iMarks fc Kroupal (|2012D . Hence, 
we assume that few-body relaxation played an important role during the initial phase of evolution 
of the ONC. In particular, three body interactions among OB stars likely led to their ejection from 
the cluster and, at the same time, to the formation of a massive object via 'runaway' physical stellar 
collisions. The resulting depletion of the high mass end of the stellar mass function in the cluster 
is one of the important points where our models fit the observational data. We speculate that the 
runaway-mass star may have collapsed directly into a massive black hole (M, > lOOM©). Such a 
dark object could explain the large velocity dispersion of the four Trapezium stars observed in the 
ONC core. We further show that the putative massive black hole is likely to be a member of a binary 
system with w 70% probability. In such a case, it could be detected either due to short periods of 
enhanced accretion of stellar winds from the secondary star during pericentre passages, or through a 
measurement of the motion of the secondary whose velocity would exceed lOkms -1 along the whole 
orbit. 

Subject headings: black hole physics — stars: kinematics and dynamics — stars: massive 



1. INTRODUCTION 

The Orion Nebula Cluster (ONC, M42) is a dense star 
cluster which is part of a c omplex star forming region at a 
distance of about 4 00pc ((Je ffries 20Q3; iSandstrom et al.l 
I2007t iMenten et all 120071 ) . Due to its relative proxim- 
ity the ONC is one of the best observationally studied 
star clusters. Its age is estimated to be < 3Myr, the 
resolved stellar mass is M c « LSOOM^ and it has a com- 
pact core of radius < 0.5pc ((Hillenbrand fc Hartmann 
I1998D . Based on the data presented by iHuff fc Stahler 
(|2006[ ) , we estimate a half-mass radius of the cluster of 
rh ~ 0.8pc. Hence, the ONC is considered to be a proto- 
type of a dense young star cluster and it naturally serves 
as a test bed for theo retical models of var ious astrophys- 
ical processes (e.g. iKroupa et al.l 120011 considered the 
dynamical evolution of ONC- type star clusters with gas 
expulsion; lOlczak et al. 2008! studied the destruction of 
protoplanetary discs due to close stellar encounters). 

Nevertheless, the morphological and dynamical state 
of the ONC is still not fully understood. Some of it's 
characteristics indicate that it has undergone a period 
of violent evolution and that the current state does not 
represent a true picture of a newl y born star cluster . 
There is a lack of gas in the cluster (iWilson et al.l[l997l ) 
which has been expelled due to the radiation pressure 
of the OB stars that reside in the cluster core. The 
effect of gas removal is likely to have led to the clus- 
ter expanding by a factor > 3 for a star formation ef- 
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ficiency < 50% (e.g. iBaumgardt fc Kroupal l2007j ) . The 
assumption that the ONC has been more compact in 
the past is in accord with the fact that this star clus- 
ter is missing wide binaries with separations > 10 3 AU 
which may have been disrupted through close thr ee-body 
inter a ctions during the i nitial compact stage (Kroupal 
120001 IParker et al.l 120091 : IMarks fc Krcmpal 12012ft . The 
process of dynamical ejections has been suggested by 
iPflamm-Altenburg fc Kroupal ()2006[ ) as a reason for the 
depletion of the cluster mass function at the high-mass 
end. This hypothesis aims to explain the observ ed deficit 
of massive stars in the ONC ([Hillenbr and 1997j) with re- 
spect to the standard IKroupa (120011) initial mass func- 
tion. Finally, observations ([Zapata et. al. 2009) have 
brought evidence that stellar disruptions occurred in the 
ONC recently. These may be a consequence of physical 
stellar collisions in the dense cluster core. 

In this paper we address the history of the ONC over 
its lifetime. We concentrate on initialy very compact 
star clusters and show that they can evolve within a few 
millions years into a state compatible with the current 
observations. In particular, we address the fact that the 
central system of OB stars, the so-called Trapezium, is 
supervirial and, at the same time, that the ONC hosts 
unexpectedly few OB stars. 

2. MODEL 

We restricted ourselves to the stage of cluster evolu- 
tion when the stars have already formed as individual 
entities which can be characterised by a constant mass 
and radius. Stellar dynamics is then driven mainly by 
gravitation al interactions . We used the numerical code 
NBODY6 (jAarsethl 120031) which is a suitable tool for 
modeling self-gravitating stellar systems with a consid- 
erable amount of binaries. 
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At the initial stage of its evolution, a substantial con- 
tribution to the gravitational field of the cluster is due 
to gas. We incorporated external gas into our model 
by means of a special type of low mass particles (m g < 
0.4 Mq) whose gravitational interaction with the rest of 
the cluster was treated directly via the N-body scheme. 
Beside that, we modified the original NBODY6 code, in- 
cluding an option for a repulsive force (with the direction 
radial from the cluster centre) acting upon the gas par- 
ticles in order to mimic the radiation pressure from the 
stars starting at a given time T ex (see Sec. 13 . 2f) . Alterna- 
tively, we also modelled the gas expulsion by an instan- 
taneous removal of the gas particles at a given t ime. 

Stel lar masses were calculated according to the lKroupal 
(|200lD mass function with an upper mass limit of 80 Mq 
(however, the mass of the most massive star is typically 
« 63 Mq). For several numerical reasons (e.g. quadratic 
growth of CPU time with the number of particles, higher 
probability of numerical errors for extreme mass ratios) 
we replaced low-mass stars (M* < m m i n ) by stars with 
mass m m i n , keeping the total mass unchanged. We have 
verified (cf. models #4 and #6 introduced below which 
have m m in = 0.5 M Q and 0.2 Mq, respectively) that our 
results are not affected considerably by this approxima- 
tion. 

The stars were cons idered to h ave finite radii i?* = 
-R (M*/M O ) - 8 (e.g. lLand 119801: this simple formula 
also fits well the radii o f zero age main seque nce stars for 
M* > M Q according to lEggleton et al.lll989l : our overes- 
timation of stellar radii for M+ < Mq does not affect our 
results, as the low-mass stars contribute only marginally 
to the merging tree) . If the separation of two stars got 
below the sum of their radii, they were merged into a 
single star. We also switched off the option of stellar evo- 
lution in the numerical code. Loss of realism is assumed 
to be negligible as we followed the cluster evolution for a 
period of only 2.5Myr. The stars in our models are not 
assigned a spectral type. For the sake of brevity, we refer 
to all stars with M* > 5 Mq as 'OB stars'. 

2.1. The 'canonical' model 

The canonical (best-fit) model of the ONC has an ini- 
tial mass of 5400 Mq, half of which is in the form of 
stars, while the othe r half accounts for gas. We consid- 
ered a Kroupal (|2001h initial mass function which, for the 
given cluster star mass, p redicts ss 50 OB stars to have 
been formed in the ONC (jPflamm-Altenburg fc Kroupal 
I2006D . In order to avoid additional statistical noise in the 
evolution of the number of OB stars, we used identical 
samples of stellar masses for all numerical realisations of 
the model. The initial half-mass radius of the cluster is 
rh 0.1 lpc, which is compatible with the initial or birth 
radius-mass relation inferred fo r star clusters ranging i n 
mass up to globular clusters (jMarks fc Kroupal 12012b . 
Positions and velocities are generated in a mass segre- 
;ated state according to the prescription of iSubr et ail 
20081) with mass segregation index S = 0.4. The al- 
gorithm by definition places massive stars in the cluster 
core and, therefore, the initial half-mass radius of the OB 
stars is w 0.05pc. In order not to place the light gas parti- 
cles at the cluster outskirts, we generated initial positions 
for stellar and gas particles separately, i.e. we had a con- 
stant gas to star mass ratio throughout the whole cluster 
at T = 0. All OB stars were set to be members of primor- 




Fig. 1. — Half-mass radius of the star cluster of our canonical 
model of the ONC (upper panel). Expansion is accelerated at 
T 0.5Myr due to gas expulsion. The total mass of the stars 
within a 3pc radius (lower panel) monotonically decreases as stars 
escape out of the observer's field of view. This curve reflects gas 
expulsion with a time delay of f» 0.5Myr which is the time the 
stars need to reach the 3pc boundary. Shaded area represents the 
la variance of the individual realisations of the model. 

dial binaries with a secondary mass M s > IMq with the 
pairing algorithm being biased towards assigning a mas- 
sive secondary to a massive primary. More specifically, 
the algorithm first sorts the stars from the most massive 
to the lighest one. The most massive star from the set 
is taken as the primary. The secondary star index, id, in 
the ordered set is generated as a random number with the 
probability density oc id~^ and max(id) corresponding to 
a certain mass limit, A/ Si , mn . The two stars are removed 
from the set and the whole procedure is repeated until 
stars with M* > M Pjm i n remain. Typically, we used as 
the minimal mass of the primary Aip iimn = 5Mq, as the 
minimal mass of the secondary M Sim j n = 1M,-., and the 
pairing algorithm index j3 = 40. For the assigned bina- 
ries we used a semi-major axis distribution according to 
Opik's law n(a) oc a -1 (jKobulnickv fc Fryer] I2007D and 
a thermal distribution of eccentricities, n(e) oc e. Ex- 
cept for the possibility of being a secondary member of 
a binary containing an OB star, we did not generate bi- 
naries of low mass stars. This limitation comes from the 
requirement of numerical effectivity and it is not likely 
to affect our results substantially. 

Let us note that the canonical model presented above 
fits into the range of t he expanding c lass of models of 
the ONC discussed by Kroupal (|2000D . Other possible 
initi al conditions would be fractal models that collapse 
(e.g. I Allison et al. I [20091 ). These require star formation 
to be synchronised across the pre-cluster cloud core to 
much shorter than the dynamical time though. 
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2.2. Dynamical evolution 

Most of the features of the evolution of a star cluster 
with strong few-body relaxation and gas expulsion can 
be demonstrated with the canonical model. In order to 
distinguish systematic effects from rather large fluctua- 
tions of cluster parameters, we present quantities aver- 
aged over 100 realisations of the model. The upper panel 
of Fig. [1] shows the temporal evolution of the half-mass 
radius of the stellar component, r^. During the initial 
phase, the cluster expands due to two-body relaxation. 
Gas expulsion that starts at T GX = 0.5Myr removes all gas 
from the cluster within a few hundred thousand years. 
From the point of view of the stars, this event leads to an 
abrupt decrease of their potential energy. Consequently, 
they can reach larger distances from the cluster centre, 
which manifests itself as an accelerated growth of at 
0.5Myr < T < lMyr. Initially weakly bound stars be- 
came unbound due to the gas expulsion and they escape 
from the cluster. Hence, reduction of the cluster mass 
also accelerates - see the lower panel of Fig. Q] where we 
plot the mass of the star cluster. We define M c as the 
sum of the mass of the stars within a sphere of radius 
Him from the cluster centre. We set r\i m = 3 pc in order 
to match typical observations of the ONC, which usually 
count stars within a projected distance of 3pc from the 
Trapezium. Note that the decrease of M c is accelerated 
with w 0.5Myr delay after the time of gas expulsion. 
This is due to the unbound stars taking some time to 
reach r\ im . Besides the escape of stars due to the gas 
expulsion, there is also continuous stellar mass loss due 
to the few-body relaxation which is capable to accelerate 
stars above the escape velocity. From T k, lMyr onwards 
the cluster expansion is again dom inated by relaxation al 
processes driving a revirialisation (jKroupa et al.ll200ll ). 

Our model indicates that the initial mass of the ONC 
must have been at least 20% larger than it is now, af- 
ter a few million years of dynamical evolution. We also 
see that, mainly due to the effect of gas expulsion, the 
ONC must have been more compact by a factor > 5 at 
the time of its birth than it is now. Hence, the initial 
relaxation time was much shorter than what we would 
infer from present-day observations. Consequently, close 
few-body interactions must have been more frequent in 
the past. In particular, we suggest that scattering of 
single stars on binaries is a process that has played an 
important role in the evolution of the ONC. This kind of 
inter action has been stud ied intensively in the literature 
(e.g. iHeggie fc Hutl [20031 and references therein). A rich 
variety of outcomes can be obtained, depending on the 
initial conditions. Nevertheless, some general results can 
be formulated. In particular, if the binary is hard, i.e. its 
orbital velocity is larger than the impact velocity of the 
third star, it typically shrinks, losing its energy which 
is transferred to the accelerated single star. Therefore, 
we expect two processes to happen in correlation: (i) 
ejections of high velocity stars and (ii) physical collisions 
of massive binaries which lead to the formation of more 
massive objects. 

The probability of collisions increases with the mass of 
the interacting stars. Therefore, numerical models by dif- 
ferent authors often show a 'runaway' process when most 
of the collisions involve t he most massive object whic h 
grows continuously (e.g. iPortegies-Zwart et alJ 12004 ). 
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Fig. 2. — Top: Mean mass of the most massive (collisional) star 
in the model regardless of its location (solid line) and mean mass of 
the most massive star within 3pc from the cluster centre (dashed 
line). Bottom: number of OB stars (solid line) and number of OB 
systems (dashed line) within 3pc from the cluster centre. 

The same is true also for our model of the ONC: The 
upper panel of Fig. [2] shows the growth of the mass of 
the most massive object, M,, due to stellar collisions; the 
number of OB stars left in the cluster is plotted in the 
lower panel. Approximately one third of the missing OB 
stars have disappeared due to merging, while the other 
two thirds have been ejected from the cluster with ve- 
locities > 10km s" 1 . The merging tree varies significanly 
among individual realisations of the particular model. In 
general, more than one merging star is formed during the 
initial m 0.5Myr. Among other processes, this is due to 
the merging of several primordial OB binaries. At later 
stages, usually one runaway-mass object dominates the 
merging process. 

In the Appendix we provide an approximate descrip- 
tion of the rate of stellar ejections due to the scattering 
of single stars on massive binaries. It gives an estimate of 
the decay time of the number of the OB stars of r > 5Myr 
for the canonical model. If we assume that the decay rate 
is linearly proportional to the number of OB stars, we ex- 
pect roughly an exponential decay of iVoB, i-e. it should 
reach half of its initial value at ~ 3.5Myr, which is in 
good agreement with the numerical results. 

A mean number of OB stars of « 28 remain at T = 
2.5Myr. This considerably exceeds the nu mber of OB 
system s observed in the ONC — the study of lHillenbrandl 
(|1997D reports only 10 stars heavier than 5M . How- 
ever, the real number of massive stars may be somewhat 
larger due to their 'hiding' in OB binary systems. Our 
models always end up with several OB stars having an 
OB companion and, therefore, a better agreement with 
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Fig. 3. — Mass of the most massive object in the cluster vs. 
the number of OB stars within a radius of 3pc. Individual points 
represent states of different realisations of the canonical model at 
T = 2.5Myr. Mass of the most massive object in the whole set 
is plotted in the upper panel, while in the lower one, only objects 
found within 3pc are considered. 



the observations is achieved if we count the OB systems 
as indicated with the dashed line in Fig. [5] Moreover, 
due to the stochastic nature of dynamical cluster evo- 
lution, different realisations of the model end up with 
quite different values of Nob and M». The final state 
(T = 2.5Myr) of all realisations of the canonical model 
is shown in Fig. [3J Although the distribution of points 
in the Nob~M, space is rather noisy, we can deduce an 
anti-corelation between these two quantities (their linear 
correlation coefficient is —0.42). In particular, all reali- 
sations that end up with Aob < 20 lead to the formation 
of a runaway-mass star of mass M. > lOOM©. In other 
words, the underabundance of the high-mass stars in the 
ONC not only indicates a period of prominent two-body 
relaxation in the past, but also the merging formation 
of a massive object which may represent an important 
footprint of the cluster's history. 

Interestingly, the term 'runaway' has a secondary 
meaning in some cases: Despite of its high mass, in sev- 
eral realisations the merging object has been ejected out 
of the cluster with velocity exceeding lOkrns -1 . These 
cases can be identified due to differences of the upper and 
lower panels of Fig. [31 The escape of the most massive 
body beyond the 3pc boundary is also responsible for the 
drops of the dashed line in the plot of the temporal evo- 
lution of AI m in Fig. [2] which shows the mean mass of the 
most massive object within 3pc from the cluster centre. 

3. DISCUSSION 



TABLE 1 

Initial parameters and final states of several models. 



id 


S 


T ox /Myr 


e 


r^/pc 


Mc /A1q 


Nob 


M,/M® 


1 


0.00 


0.5 


0.50c 


0.56 


2183 


35.1 


102.9 


2 


0.00 


0.7 


0.50c 


0.48 


2061 


28.1 


138.1 


3 


0.25 


0.7 


0.50c 


0.51 


1980 


26.2 


140.9 


4 


0.40 


0.5 


0.50c 


0.59 


1895 


28.0 


140.1 


5 


0.40 


0.5 


0.50c 


0.60 


1938 


27.3 


138.5 


6 


0.40 


0.5 


0.50c 


0.63 


1887 


28.3 


144.2 


7 


0.40 


0.5 


0.50c 


0.66 


1660 


26.0 


142.7 


8 


0.40 


0.5 


0.50v 


0.47 


2197 


29.7 


122.1 


9 


0.40 


0.7 


0.50c 


0.62 


1895 


26.5 


152.1 


10 


0.40 


0.7 


0.50c 


0.58 


1882 


27.0 


144.2 


11 


0.40 


0.7 


0.50c 


0.67 


1636 


26.1 


132.6 


12 


0.00 


0.7 


0.33c 


0.52 


1717 


24.2 


173.7 


13 


0.25 


0.5 


0.33v 


0.51 


2109 


32.1 


124.2 


14 


0.40 


0.5 


0.33v 


0.50 


2196 


33.1 


112.1 


15 


0.40 


0.5 


0.33c 


0.77 


1300 


23.2 


171.0 



Common initial parameters of all models are the initial mass of 
the stellar component, M c = 2700Mq which implies Nob(T = 
0) = 50 and M.(T = 0) « 63M Q . The initial half-mass radius 
varies slightly throughout the models, but is generally a O.lpc. 
The gas particles are of mass 0.2Mq for models 6 and 10; in all 
other cases m gas = OAMq is considered. In models 7 and 11, gas 
particles were removed instantaneously at T = T ex ; in all other 
cases, T ex is the time when the repulsive force was switched on. 
S is the mass segregation index as defined in I Subr et al.l H2008T ) . 
e = M C /(M C + Mgas) is the star formation efficiency; suffix 'c' 
stands for constant e throughout the whole cluster, while 'v' means 
variable e (decreasing outwards). Opik's distribution of semi- major 
axes of the primordial binaries is considered in all models except 
for #5, where we set n(a) = const. The canonical model has id 6. 

The final state of the canonical model presented in the 
previous section matches basic observables of the ONC 
quite well, in particular its mass and half-mass radius. 
The mean number of OB stars is somewhat larger than 
what is observed in the ONC, nevertheless, some of the 
realisations reach Aob < 20. Hence, we consider this 
model to be a realistic representation of the ONC. We 
have investigated several tens of different models with 
different values of the parameters including the index 
of the initial mass segregation, initial half-mass radius, 
mass of the cluster, star to gas mass fraction and the 
time of the gas expulsion. Models whose final state can 
be considered compatible with the current state of the 
ONC, at least in some of their characteristics, are listed 
in Table [T] 

3.1. Nob — M t anticorelation 

As mentioned above, the canonical model indicates an 
anticorelation between the final number of the OB stars 
and the final mass of the runaway-mass star. In Fig. @]we 
plot the mean mass of the most massive object vs. the 
mean number of remaining OB stars for all models listed 
in Table [1] Now, each point represents an average over 
several tens of realisations of the particular model and a 
Nob — M, anticorelation becomes evident. We attribute 
this relation to the fact that both mechanisms of OB star 
removal, i.e. physical collisions and ejections, are driven 
by a common underlyig process of close three-body in- 
teractions (see Appendix). Naturally, their importance 
grows with increasing stellar density. Hence, the models 
that stay more compact during the course of their evolu- 
tion are located in the top-left corner of the graph. They 
better fit the observations from the point of view of the 
number of OB stars, however, either their half-mass ra- 
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Fig. 4. — Mean mass of the most massive object in the cluster vs. 
the mean number of OB stars (both within 3 pc) at T = 2.5Myr for 
several different models of compact star clusters listed in Table [T] 
The canonical model is represented with the X sign; full circle 
stands for the initial state which is common for all models. 

dius or total stellar mass at T = 2.5Myr is too small, i.e. 
not consistent with the observations. 

3.2. Gas expulsion 

In most of our models, the gas particles were blown out 
of the cluster due to a repulsive external force oc r/r 3 
centered on the cluster core. The strength of the force 
was set such that the gas particles were accelerated to ve- 
locities » 8km s -1 on the time-scale of w 0.1 Myr. This 
approach is the most realistic one. keeping the continu- 
ity equation fulfilled, but, at the same time, it is nu- 
merically the most expensive. The alternative method 
with instantaneous removal of the gas particles at T ex 
leads to similar results. The clusters integrated with this 
method typically result in somewhat lower total mass and 
larger half-mass radius at the final time. Let us note, for 
completness, that in the literature, another method that 
mimics gas remov al in numerical models o f star clusters 
is also used (e.g. iGever fc Burkertl 120011 : iKroupa et al.l 
I2001L iBaumgardt fc Klessenl 120111) . It is based on the 
time-variable smooth external potential that represents 
the weakening gravitational potential of the gas. It has 
been found that the two approaches (i.e. modelling gas 
via particles vs. ext ernal potential) lead to nearly indis- 
tinguishable results (|Gever fc Burkertll200lD . In order to 
check for a possible bias of our gas particles approach 
which may stem from two-body relaxational processes, 
we have integrated two additional models with very low 
masses for the gas particles. Due to the large num- 
ber of particles in these integrations, we have followed 
their evolution only up to the time of the gas expulsion, 
T cx = 0.5 Myr. Outcomes of these models in terms of the 
mean number of remaining OB stars, Nob, and the mass 
of the runaway-mass star, M # , are presented in Table [2] 
Besides just giving the mean values of -/Vob and M., we 
have also performed a Kolmogorv-Smirnov test compar- 
ing the sets of Nob and M. coming from these models 
with the canonical one. As can be seen, there is no ap- 
parent dependence of the results on m gas ranging from 
0.03 to 0.4 Mq. Hence, we conclude that the two-body 
relaxation due to the gas partiles does notsignificantly 
affect our results. 

The evolution of the radiation pressure in the real ONC 
has been definitely more complicated. According to ob- 



TABLE 2 

Comparison of models with different masses of gas 
particles 



m g 


Nob 


M m 


Ar run 


P KS,Nor 


-Pks.m. 


comment 


0.40 


39.7 


125.1 


10 


95% 


42% 


model #4 


0.20 


39.8 


124.9 


100 






canonical model 


0.05 


38.2 


128.5 


14 


82% 


87% 




0.03 


39.9 


126.7 


11 


32% 


75% 





Mean values of the number of OB stars within the radius of 3 pc 
from the cluster core, Nqb, an d mass of the runaway-mass star, 
M m , for models with different values of m g . All other parameters 
of the models, including the total mass of gas, are identical to the 
canonical model. Pks,* 1s the Kolmogorov-Smirnov probability 
that the null hypothesis (assuming the particular data set comes 
from the same distribution as the canonical one) is valid. Note 
that the value of Pks,* ~ 5% is usually considered to be a limit 
below which the null hypothesis should be rejected. 

servations (jHuff fc Stahlerl f2"006). the star formation in 
the ONC was continuous with the most massive stars 
starting to be formed no more than 2 Myr ago. Setting 
in our models to either 0.5 or 0.7 Myr appears to have 
only a marginal impact on the final state of the cluster. 
Hence, we assume that a more complicated temporal pre- 
scription of the gas removal would not affect our results 
considerably. 

3.3. The runaway-mass star or black hole 

Under the assumption of an invariant canonical IMF 
and a small initial cluster radius, the formation of the 
runaway-mass object via stellar collisions appears to 
be an inevitable process that, together with the high- 
velocity ejections, decreases the number of OB stars in 
a dense star cluster. It's possible detection would defi- 
nitely strongly support the scenario of the ONC history 
presented in this paper. 

The current state of the runaway-mass star is, how- 
ever, not clear. The massive star would definitely have 
undergone a fast internal evolution. While successive 
merging events may lead to very fast mass growth, stellar 
winds act in the opposite way. The state of the runaway- 
mass object after « 2 Myr of evolution depends on the 
(dis)balance of these two processes. In the literature, 
quite different predict ions on this subject can be found. 
iGlebbeek et al.l (|2009D suggest that winds should work 
in a self-regulatory manner, limiting the maximum mass 
and, consequently the star's lifetime. Oth er authors (e.g. 
ISuzuki et alll2007l : lPauldrach et al.H2011[ ) state that stel- 
lar winds of the runaway-mass star will not be able to 
terminate its growth and it may then collapse to a mas- 
sive black hole. There is no observational evidence for a 
star with a mass above 50M Q in the ONC which could be 
interpreted as a runaway-mass object. The most massive 
stellar member of the ONC, 1 C, is a binar y with a total 
mass of about 45A/ Q and mass ratio sw 0.25 (jKraus et al.l 
l2009t) , i.e. the primary has mass < 35 Mq. It does not 
exhibit a stellar wind strong enough to reduce its mass 
by several tens of Mq. Hence, our scenario assumes a 
low efficiency of stellar winds, i.e. continuous growth of 
the runaway-mass star which forms a massive black hole 
at the end of its lifetime. There is no general consensus 
whether the collapse of a star more massive than 150Af Q 
will be followed by an ejection of most of its mass. As 
there is no evidence for a supernova remnant in the ONC, 
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Fig. 5. — Histogram of the semi-major axes (a, empty boxes) 
and pericentre distances (r p , full boxes) of the binaries involving 
the runaway-mass star in the canonical model. A total of 100 runs 
were examined, out of which 70 harbour the runaway-mass star 
with a stellar companion (three of them with a > 1000AU). 

we assume most of the mass of the runaway-mass star, if 
it is present, to have collapsed directly into a black hole. 

The putative black hole in the ONC would be de- 
tectable only through an interaction with its environ- 
ment. One possibility is an accretion of surrounding gas 
which could be detectable through X-ray radiation. The 
main source of the gas in the Trapezium region are stel- 
lar winds from the remaining OB stars. We estimate 
that they produce at most lO^Moyr -1 within the cen- 
tral 0.2pc. However, the black hole is likely to capture 
only the gas that falls within its Bondi radius « 100 AU, 
i.e. the accretion rate could be < 10~ 15 MQvr . Such 
a highly sub-Eddington accretion will not be detectable 
in the region confused with several X-ray young stellar 
sources. 

A better chance for detecting the black hole would be 
given if it is a member of a binary. This indeed appears 
to be the case for about two thirds of the realisations of 
the canonical model. However, the typical separation of 
the black hole and the secondary appears to be from sev- 
eral tens to hundreds of AU (see Fig. [5) • Even with the 
relatively large eccentricites achieved, in no case does the 
secondary star fill its Roche radius at the pericentre of its 
orbit, i.e. we do not expect the black hole to be a mem- 
ber of a mass transferring system. Still there would be a 
chance for a detectable accretion provided the secondary 
star is massive (M+ > 10 Mq). In such a case, its stel- 
lar wind may produce short periods of enhanced activity 
during pericentre passages, but the corresponding lumi- 
nosity is rather uncertain. Let us consider the well kown 
Cyg-Xl system as a kind of template. Having a black 
hole mass of « 10 Mq and separation from the secondary 
star of sa 0.2 AU, i t has a luminosity w 2 x 10 37 ergs _1 
([Wilms et al.l 120061 ) which corresponds to « 10~ 4 of the 
rest-mass energy of the winds emitted by the secondary 
star (« 10 -6 Mq yr _1 ). Assuming the same effectivity in 
conversion of the stellar winds into radiation, we estimate 
the luminosity to be similar to the Cyg-Xl system, i.e. 
« 10 4 Lq, for separations of the order of 1 AU. Consid- 
ering the density of the stellar wind to decrease with the 
sqaure of the distance, we obtain a rough estimate of the 
peak luminosity L « 10 4 (r p /l AU) _2 L Q . Acording to 
the mass of the putative black hole, the maximum of the 
emitted radiation is expected to lie in the X-ray band. 
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Fig. 6. — Velocity dispersion a within a specified radius r at the 
final state of the canonical model. For distinguished bound multi- 
ple systems, the centre-of-mass velocity is counted instead of the 
proper velocities of the individual components. Solid line repre- 
sents the average over all realisations; thin dashed line corresponds 
to half of the computations yielding more massive (M. > 120 Mq) 
runaway-mass object, while the dotted line stands for the com- 
putaions with lower M m ■ 

The accretion events should be recurrent with a period 
ranging from years to hundreds of years. A shorter pe- 
riod generally implies a smaller r p , i.e. larger values of 
the peak luminosity. 

3.4. Velocity dispersion 

Another way of an indirect detection of the massive 
black hole in the ONC lies in velocity dispersion measure- 
ments. It has been stated by several authors that the core 
of the ONC (r < 0.25 pc) is dynamically hot with a ve- 
locity dispersion > 4kms~ 1 (e.g. iJones fc Walkerlll988l : 
Ivan Altena et al.l KL988I : iTobin et all 12009( 1. Themissing 
stellar mass required for virial equilibrum in the inner- 
most region of the ONC was estimated to be > 2000 Mq 
which exceeds th e observed mass sa 200 Mq by an or- 
der of magnitude. Kroupa et al.l (|200lH demonstrate that 
the globally super- virial velocity dispersion is readily ob- 
tained if the ONC is expanding now after expulsion of its 
residual gas, but some of the missing mass in the inner 
region could be attributed to the invisible remnant of the 
runaway-mass star. 

As we show in Fig. [5] our canonical model gives a ve- 
locity dispersion > 4km s" 1 in the inner 0.25pc, which 
is in accord with the observational data. In order to 
inspect the influence of the mass of the runaway-mass 
object on the velocity dispersion, we have divided the 
100 realisations of the canonical model into two groups 
of 50 according to the mass of the most massive object 
that remained in the cluster. It appears that the clus- 
ters with M, > 12OM have a velocity dispersion within 
0.25pc somewhat larger than those with M, < 120Af Q , 
nevertheless, even within the latter group we have a > 
3.5km s~ , i.e. an apparently super- virial core region. 
Together with the fact that none of our numerical ex- 
periments led to M, > 2000 Mq, virtually required by 
the condition of virial equilibrium, this indicates that 
the large velocity dispersion can be only partly due to 
the hidden mass of the runaway-mass object. Detailed 
analysis of our models shows that there are two other 
reasons for having large velocity dispersion in the core. 
First, according to our models, the ONC is a post-core- 
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TABLE 3 



Trapezium 


STARS DATA 


name 


M/Mq 






18.9 


33.4 ± 2 


e x s 


7.2 


24.0 ± 2 


e 1 ^ 


44 ±7 


23.6 


e 1 D 


16.6 


32.4 ± 1 



Stella r masses of Q 1 ^!, S 1 ^ and @ 1 D are taken from IHillenbrancD 
Jl997f> . Radial velocity data of G 1 A and Q^B come from the C.D.S. 
- SIMBAD database; data of S 1 D follows IVitrichenkd (2002T). 
Both m ass and radial velocity of S 1 C are according to lKraus et al.l 
(EH). 

collapse cluster with centrally peaked density. This leads 
to a deep gravitational potential well which allows the 
central velocity dispersion to be larger than what would 
be expected e.g. for the Plummer profile. Second, in 
our models and perhaps also in the real ONC, there are 
likely to be present several not identified wide binaries 
with orbital velocities exceeding lOkms -1 which affect 
the velocity dispersion. 

Fig. [5] indicates that the presence of the dark mas- 
sive body influences remarkably (i.e. by more than 
1 km s~ ) the velocity dispersion in a small region of a 
radius of rs 0.1 pc. This is in accord with an analyt- 
ical estimate of the influence radius of the black hole, 
r h = GM./a 2 , which for M. = 100 M Q and a = 3kms _1 
gives rh ~ 0.05 pc. Due to the small number count of 
stars in such a small region, standard statistical meth- 
ods do not represent a robust tool for determining the 
presence of a dark massive body. Nevertheless, hav- 
ing this caution in mind, let us discuss the velocity 
dispersion of the compact core of the ONC, the four 
Trapezium stars. According to the published observa- 
tional data (see Table [3} , the line-of-sight velocity dis- 
persion of the Trapezium is o\ a& « 4.6kms~ which im- 
plies a 3D velocity dispersion a 7.9 km s -1 . Consid- 
ering a sphere of radius rx ~ 0.025 pc covering the four 
Trapezium stars, we obtain an estimate of the enclosed 
mass required for the system to be in virial equilibrium: 
A^bind ~ rT<j 2 /G « 350 Mq. Taking into account the 
stellar mass of the Trapezium (w 90 Mq), we find it to be 
an apparently dynamically very hot system. Our models 
suggest that the virial equilibrium state can be achieved 
by the presence of the runaway-mass object. 

Yet another, and statistically better founded, argu- 
ment for a hidden massive body comes from compari- 
son of the kinematical state of the Trapezium and sim- 
ilar subsystems found in our numerial models. We 
have implemented an algo rithm similar to that used by 
lAllison fc Goodwml (|2011D for detection of Trapezium- 
like systems. In particular, for each OB star we de- 
termined its three nearest OB neighbours (either sin- 
gle OB stars, or bound systems containing at least one 
OB star). The set of stars was considered Trapezium- 
like if all members had the same OB neighbours. In 
order to obtain a statistically significant sample, we ex- 
amined the last 20 snapshots (covering the period from 
« 1.5 Myr to 2.5 Myr) of 100 independet realisations 
of the canonical model. We found Trapezium-like sys- 
tems in a third of cluster snapshots. The mean value 
of the velocity dispersion of these systems was found 



to be 7.5 kms -1 . The runaway-mass star was ex- 
cluded from the search algorithm. Furthermore, we 
have distinguished systems according to their centre of 
mass distance to the runaway-mass object. In approx- 
imately two thirds of the Trapezium-like systems, the 
runaway-mass object was found within a distance smaller 
than max( Ar^ ) , with Ary denoting separations of indi- 
vidual members. The mean velocity dispersion of the 
Trapezium-like system of this subset was m 9.0 kms - , 
while it was only w 4.3 km s -1 for the remainig cases. 

Finally, let us remark that kinematical evidence of a 
massive black hole in the ONC may also lie in the w 70% 
probability that it is a member of a binary (see Sec. I3.3|) . 
The typical velocity of the secondary star should be 
> lOkms -1 , i.e. considerably exceeding the observed 
central velocity dispersion. 

4. CONCLUSIONS 

We have carried out extensive modelling of the dy- 
namical evolution of compact young star clusters, aimig 
to reconstruct the history of the ONC. Assuming that 
the ONC underwent primordial gas expulsion, we have 
shown that the ONC must have been sev eral time s more 
co mpact than it is now in agreement with lKroupal (l2000h 
and lKroupa et all (|2001l ). This implies that its two-body 
relaxation time was considerably shorter in the past and, 
consequently, that close few-body interactions between 
massive stars were rather frequent. We have concen- 
trated on two significant tracks of such interactions that 
lead either to high-velocity ejections of massive stars 
from the cluster or to their physical collisions that are 
likely to lead to the formation of a massive 'runaway' 
merging object. Both of these processes decrease the 
number of massive stars in the cluster. Hence, the ob- 
served significant lack of massive OB stars in the ONC 
further supports our assumption that this star cluster 
has undergone a relaxation dominated period since its 
birth. 

We have shown that it is not unfeasible that the centre 
of the ONC harbours a black hole of mass > 100M© as 
the remnant of the massive runaway-mass star. Its pres- 
ence could be revealed either by episodic accretion events 
or by a thorough kinematical study of the innermost 
« 0.05 pc region of the ONC, where the hidden mass 
would increase the velocity dispersion above 5km s _1 . In 
particular, we have shown that the observed velocity dis- 
persion of the Trapezium system can be achieved by the 
presence of an object of mass « 150 Mq which is fully 
consistent with the hypothesis of the runaway-mass ob- 
ject. Our model also shows that the apparent super- 
viriality of the central w 0.25pc can be explained as a 
natural attribute of a post-core-collapsc dynamical state 
of a star cluster with a considerable binary fraction. 

The possible detection of the remnant of the merger 
object could bring new light into several fields of con- 
temporary astrophysics. First, it would confirm the hy- 
pothesis that star clusters similar to the ONC are being 
formed very compact with a ha lf-mass radius of the or - 
der of a few tenths of a parsec (jMarks fc Kroupall20ll . 
Second, and probably more importantly, it would have 
important implications for the evolution of very massive 
stars and merger products for which we have only a lim- 
ited understanding now. 
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APPENDIX 
THEORETICAL PREDICTIONS 

Unlike two-body relaxation, interactions among three or more stars allow considerable energy transfer from one 
component to another. Therefore, multiple-body scattering is essential for both processes (merging and ejections) that 
lead to the reduction of the number of massive stars in the cluster. In spite of the chaotic nature of the dynamical 
evolution of the star cluster core, it is possible to derive a raw esti mate of the rate of OB stars ejections and collisions. 

Our approximation of three-body scattering follows the work of lPerets fc Subrl ()2012D : Consider an interaction of a 
binary of mass Mb = Mi + M 2 and a single star M*. Large accelerations of the single star are assumed to occur when 
it passes around one of the binary components within the semi-major axis, a, of the binary. The cross section of the 
interaction is assumed to be determined by gravitational focusing: 
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S(a) 



2nGM B a 



(Al) 



where v c is the characteristic stellar velocity in the cluster and G stands for the gravitational constant. 

The energy transfer estimate is based on the approximation that the impacting star moves along a hyperbolic orbit 
around Mi perturbed by M 2 . The typical perturbing force is 

GM 2 M+ 



F ; 

and it acts along the star's trajectory segment of length ps a. Hence, the energy transfer to the star is 

. _ GM 2 M+ GM B M+ 



(A2) 



(A3) 



Here we for simplicity assume the mass of the binary to be of the same order as the mass of the secondary, AIb ~ M 2 . 
In the cases when Ai?* is much larger or at least comparable to the star's energy before the interaction, it will be 
accelerated to 



2GM B 



(A4) 



Finally, let us assume the distributions of the binary mass, semi-major axis and eccentricity in a simple power-law 
form, 



n(a, e, M B 



2Aea~ 1 M^ a 



(A5) 



Catch me if you can 



9 



in the interval a G (a m ; n , a m ax); e G (0,1) and Mb G (M m ; n , M max ), i.e. n(a, e, Mb) da de c!Mb is the number of 
binaries with semi-major axis, eccentricity and mass in (a, a + da), (e, e + de) and (Mb, Mb + dM B ) respectively. The 
normalisation constant then reads: 

r^feV^fe) for a=l > 

-4 = 1 / ( (A6) 



In 



1 f 2»«0 i_ Q for a ^ 1 



ff/Gff- VELOCITY EJECTIONS 

The frequency of scattering events of a star with velocity v c on binaries with number density n B and semi-major 
axis a is 

2TTGM B n B a 

v = L n B u c = , (A7) 

where we assumed the binary cross-section to be determined by the gravitational focusing, i.e. X w 2i:GM B a/v ^. The 
mean frequency of ejections can be obtained via integration of (|A7j) weighted by the distribution function (|A5[) : 

/■Afmax /"Glim 

jy e = / dM B / da^(a,M B )n(o,e,M B ) . (A8) 

The upper limit of the semi-major axis, aii m , has to be set such that only interactions that lead to acceleration above 
the escape velocity from the cluster are considered, i.e. 



, n 2GM C , k . 

Vacc(aiim) = V csc fH i , (A9) 

V 

where r c is a characteristic radius of the cluster. Combining ()A9|) and (|A4[) gives an m « r c M B /M c . For the canonical 
model with initial M c = 5400M o and r c = 0.1 lpc and for 5M Q < M B < 100M Q we obtain a Um > 100 AU which is the 
upper limit of the semi-major axis distribution used in our numerical integrations. Hence, let us for simplicity assume 
a>\im = flmax = const, which yields 



27rGn B ^max ^min ^max -^min r ^ 

T~i / \ , / , f 7T7 r lor a = 1 

u c ln(a max /a min ) ln(M max /M min ) 

2ttGhb a max — a min 1 — a M,„ ax — M mll ^ ^ a 7^ 1 

u c ln(a max /a min ) 2 - a Mml" - M^7 n Q 



(A10) 



Considering an initial OB binary density tib ~ 10 pc 3 and velocity dispersion v c w 10km s , a G (0.1AU, 100AU) 
and Mb G (5M , 100M q ) with a Salpeter mass function (a = 2.35) we obtain D c O.lMyr -1 . Hence, the mean time 
for an OB star to undergo scattering on a massive binary that leads to its ejection from the cluster is r = \ jv e lOMyr 
for the canonical model. 

STELLAR COLLISIONS 

Unlike in the case of stellar ejections, we assume that stellar collisions are caused by binary shrinking due to successive 
three body interactions. We further assume all the massive binaries to be hard, i.e. their interaction with the third 
body leads to the acceleration of the impact star and growth of the binding energy of the binary. The frequency of 
the events of scattering of a star of mass M* on the massive binary is v 1 — Sw c n(M v )dM A ., where n(M*)dM* is the 
number of stars of given mass per unit volume. The rate of the energy transfer from the binary to the stars of mass 
in (M+, Ah + dM*) is 

dE , 2-kG 2 MIM+ , s s 

— « AEy = 3_i n(M*)dM* . All 

di v c 

The net rate of energy transfer to stars within the whole mass spectrum gives 

^ « / = 2 ^ M ^ . (A12) 

dt J v c v c 

The stars collide when the binary separation at pericentre becomes less than the sum of the stellar radii. In terms of 
binding energy and with approximation i?*(Mi) + i?*(M 2 ) w 7?*(Mb), the condition for collision reads: 

E co n « '-l^ML . (A13) 
8 iJ*(M B ) V ' 

The time required to grow the binding energy from the initial value Eq « ^GM^/clq to E co a is 



icoll ~ (Eo 



dE 



dt 



-1 



1 - e 1 



8ttG>* \R*(M B ) a 



(A14) 
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Integration of t co u over the whole parameter space of binaries with the distribution function (|A5j) gives the characteristic 
time of stellar collisions: 



T coll ' 



da 



da 



M„ 



M mi; 
M„ 



dMv 



l-R i ,(M B )/a 



de 



-«( 1 



1 



1 Ri 

1 2 



1 

a 



ra(a, e, M B ) 



3a 3 



(A15) 



where we denoted i?B = -R*(Mb); the upper limit on the eccentricity is set such that the binaries are not collisional 
initially. Assuming R B 3> a, we omit the two least significant terms in ()A15|) and perform the integration over a: 



Tcoll 



Mr, 



Mr, 



dM B Mb ' 



3i?B 



In 



(A16) 



Finally using the notation M = M/Mq we obtain for R+(M) = Rq M 



o.s 



167rGp*i?0 



12 V 



-0.8 



Mr 



-0.8 



In" 



i?,: 



Tcoll 



lGTrGjO^i?© 



1 



0.2-q 



0.2-q 



0.6 - 3a Ml 



Mi 



R, 



Rp> 



In" 



In" 



for a = 1 , 



for 



(A17) 



a ^ 1 . 

The canonical model has an initial central density p 4 «2x 10 e M Q pc -3 which gives r co ii ~ 40 Myr, i.e. we estimate 
stellar collisions to be roughly a factor of 4 less efficient process for OB star removal than three-body scattering. This 
is somewhat less than what our numerical experiment shows. This discrepancy can be due to several reasons. Most 
important is probably the fact that stellar collisions may occur sooner due to perturbations which lead to a strong 
growth of orbital eccentricity. Indeed, the numerical experiments show that most of the collisions occur with e > 0.99. 
Hence, the above given derivation can only serve as an order of magnitude estimate. 

Altogether, the three-body interactions are expected to decrease the number of OB stars in the canonical model of 
the ONC by a factor of 2 on a time-scale of « 5 Myr which is in accord with the results of the numerical experiment. 



